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Abstract 


The Phase Space is a powerful tool for representing and reasoning about the qual¬ 
itative behavior of nonlinear dynamical systems. Significant physical phenomena of 
the dynamical system - periodicity, recurrence, stability and the like - are reflected 
by outstanding geometric features of the trajectories in the phase space. Success¬ 
ful use of numerical computations to completely explore the dynamics of the phase 
space depends on the ability to,(l) interpret the numerical results, and (2) control 
the numerical experiments. This paper presents an approach for the automatic re¬ 
construction of the full dynamical behavior from the numerical results. The approach 
exploits knowledge of Dynamical Systems Theory which, for certain classes of dynam¬ 
ical systems, gives a complete classification of all the possible types of trajectories, 
and a list of bifurcation rules which govern the way trajectories can fit together in 
the phase space. These bifurcation rules are analogous to Waltz’s consistency rules 
used in labeling of line drawings. The approach is applied to an important class of 
dynamical system: the area-preserving maps, which often arise from the study of 
Hamiltonian systems. Finally, the paper describes an implemented program which 
solves the interpretation problem by using techniques from computational geometry 
and computer vision. 
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Everything should be made as simple as possible, 

BUT NOT SIMPLER. 
- A. Einstein 


1 INTRODUCTION 


Qualitative Physics is a young field. Progress is made when researchers formalize 
and implement their understanding of how certain qualitative reasoning tasks, such 
as prediction of future behavior, and explanation of how the behavior comes about, 
are being performed in particular problem domains. Two domains, among others, 
have received much attention: circuit analysis and design in the engineering domain, 
and simple boilers and fluid flow in commonsense physics. Early works in Qualitative 
Physics primarily dealt with incremental deviation from equilibrium states where time 
evolution is not explicitly considered [4]. More recent works attempt to extend DeK- 
leer’s qualitative algebra and incremental analysis to handle time-varying behavior 
[5,15,14,11]. 

A common characteristic of the examples examined so far is this: the physical 
dynamic systems studied in these disciplines can often be divided into different op¬ 
erating regions each of which can be modeled to a first approximation by simple low 
order linear differential equations, or by first order nonlinear differential equations. 
The behavior of linear systems is particularly simple: the complete input-output be¬ 
havior can be summarized in a single system transfer function. Consequently, if the 
response to one type of input is known, no more information is needed to determine 
responses for other input signals. 

The situation in a nonlinear system is completely different: essential changes in 
the qualitative behavior of the system may occur as the amplitude of the input signal 
changes, or as the starting conditions are varied. More importantly, nonlinear systems 
have a far richer spectrum of dynamical behavior. Simple equilibrium points, periodic 
and quasiperiodic motion, limit cycles, chaotic motion as unpredictable as a sequence 
of coin tosses - these are some of the behavior found in a typical nonlinear system. 

Unfortunately, these nonlinear characteristics do not show up in first order non¬ 
linear differential equations. This is because the continuity and (local) uniqueness of 
flow severely constrain the kind of behavior possible on the real line: the flow either 
tends towards an equilibrium, or goes off to in fini ty 
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In this research, I therefore propose to look at dynamical systems - those typically 
encountered in Physics - to provide a new source of examples for investigation into 
the fundamental issues of descriptive language, style of reasoning, and representation 
techniques in qualitative reasoning about nonlinear dynamical systems. Specifically, I 
will consider two-dimensional discrete dynamical systems defined by area-preserving 
maps containing a single control parameter. The study of area-preserving maps - 
transformations of the plane which preserves area - began with the venerable prob¬ 
lem of the stability of the solar system. I choose to investigate this simplest non-trivial 
type of conservative system because many important problems in physics - the re¬ 
stricted 3-body problem, orbits of particles in accelerators, and two coupled nonlinear 
oscillators, just to mention a few — can be reduced to the study of area-preserving 
maps. 

To illustrate the vocabularies that a physicist uses to describe the qualitative 
behavior of nonlinear systems, we can examine the following description of one of the 
phase portraits obtained from the results of numerical experiments with the quadratic 
map (see next section): 


Figure la, represents a number of trajectories (sequences of points) for cos a = 
0.4 ... a regular structure in a neighborhood of the elliptic fixed point at the 
origin, and farther away a chaotic zone. In many places the plotted points are 
so dense that they give the illusion of a continuous curve. Near the origin, 
the “curves” are almost circular. As we move outward, the curves become 
distorted. Just inside the outermost regular curve lies a chain of six closed 
curves, or “chain of islands”. Successive points of a trajectory jump from one 
island to another by application of the mapping. Finally, as the curves become 
more distorted, there is a sudden break-up and the set of points no longer lies 
on a curve, but seems rather to fill a two-dimensional region. 


The description is an edited form of a passage from [8, pages 99-100]. What is 
striking in this description is that the language is entirely geometric: it contains terms 
like fixed point, continuous curve, chain of islands, two-dimensional region etc. The 
phase portrait is completely characterized by the spatial relationships among these 
geometric objects. 

A goal of this research is to develop a clear understanding of the geometric lan¬ 
guage displayed in the above description. In particular, I want to understand how 
representation and reasoning techniques based on the geometric language can be used 
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as a basis for a program to perform automatic numerical experiment control and in¬ 
terpretation. 

The paper is organized as follows. Section two defines the task. Section three 
reviews the knowledge — qualitative dynamics and bifurcations — required for the 
task. Section four describes my approach to the experiment control problem. This 
part has not been implemented. Section five, the main result of this paper, shows the 
major pieces of an implemented program which solves the interpretation problem. 


2 THE TASK 


Given an one-parameter area-preserving map defining a discrete dynamical system, 
I am interested in describing the qualitative long-term, or asymptotic, behavior of 
the system over its entire operating range. Specifically, if U is an open subset of the 
phase space X C 3£ x 3£, and J C an interval over which the control parameter 
varies, I want my program to automatically generate (1) a family of phase portraits 
describing the main dynamical properties of the map for all initial conditions in U 
and parameter values in J, and (2) a summary graph which shows a cross-section of 
the map along certain axis of symmetry as a function of the parameter. 

To explore the complete dynamics of a nonlinear system over a large region of the 
phase space and parameter space is a fairly typical problem in the physics literature. 
A good illustration of this task is provided by Henon’s well-known paper, “Numerical 
Study of Quadratic Area-Preserving Mappings” [9]. The goal of Henon’s paper is to 
provide a description of the main properties of the quadratic map: 

x n+ 1 = COS a — (y n — xl) sin a 

2/„+i = # n sin a + (y n — x„) cos a (1) 

where x and y are the state variables, and a is the control parameter. The main 
results of Henon’s paper are shown in Figures l(a)-(f), which display the output of 
many numerical simulations, and Figure 2 which summarizes the dynamics of the 
map in a two-dimensional plot. Dots on the plot indicate regions of invariant curve, 
circles indicate islands, and blanks indicate chaotic zones. 

The simplest approach to this problem is the brute force method: it divides the 
phase space and parameter space into small grids and tries every possible combi¬ 
nations of initial conditions and parameter values. A simple calculation will show 
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Figure 1: A partial list of phase portraits from numerical experiments, (a) a = 1.16 
(b) a = 1.33 (c) a = 1.58 (d) a — 2.0 (e) a = 2.04 (f) a = 2.21. Dashed line: axis of 
symmetry. 

that this involves an enormous amount of computation. For instance, if we choose 
a uniform grid size of 0.01, we have to compute approximately 300 x 300 x 600 = 
54 million orbits. Assuming, on the average, 0.02 second is needed to compute a 
trajectory of 500 points, it will take over 300 hours of computation time to compute 
all the trajectories. 

The brute force method suffers from two serious problems. First, it is grossly 
inefficient because most of the phase portraits computed will be qualitatively the 
same. Second, it is not reliable because there is always the danger of missing some 
important qualitative features when the change occurs at a resolution finer than the 
grid size. 
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Figure 2: Summary picture for all values of a, and initial points on the axis of 
symmetry. Dots: regions of invariant curve. Circles: islands. Blanks: chaotic zones. 
Dashed curves: unstable invariant points. 

A physicist often does much better than this. Figure 3 represents a flow-chart of 
what an expert physicist does during the numerical experiment. The flow-chart has 
two nested loops. The outer loop involves deciding when to stop the experiment; the 
inner loop, when to move on to next parameter value. Controlling what experiment 
to do next, and interpreting the results of the simulation - these are the two most 
important decisions the experimenter has to make. 

The task of behavior prediction can now be summarized as this: to develop a 
picture of all possible solutions to the dynamical system from a limited amount of 
numerical experiments at a limited number of initial conditions and parameter values. 
The key observation is that knowledge of qualitative dynamics and their geometric 
manifestations in the phase space provides a strong constraint on the type of behavior 
possible. As we will see in the next section, this constraint translates into a dramatic 
reduction of the amount of search required to find those combinations of states and 
parameter values that lead to “interesting” phase portraits. 


3 WHAT IS THE KNOWLEDGE? 
3.1 Terminology 


The purpose of this section is to introduce some concepts and definitions from Dy¬ 
namical Systems Theory [10]. A dynamical system consists of two parts: (1) the 
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INPUT: Equfction* of mnlion 
R*nge of Initial States 



Figure 3: Flow-chart which describes the process of experimenting with a dynamical 
system. 

system state, and (2) the evolution law. The system state at any time t 0 is a mini¬ 
mum set of values of variables , x n } which, along with the input to the system 

for t > t 0 , is sufficient to determine the behavior of the system for all time t > t 0 . The 
variables which define the system state are called state variables. The conceptual 
n-dimensional space with the n state variables as basis vectors is called the phase 
space. A state vector is a set of state variables considered as a vector in the phase 
space. As the system evolves with time, the state vector "traces out a path in the 
phase space; the path is called an orbit or a trajectory. Finally, a phase portrait 
is a partition of the phase space into orbits. 

The evolution law determines how the state vector evolves with time. In a finite 
dimensional discrete time system, the evolution law is given by difference equations. 
The difference equation is specified by a function / : X —► X where X is the phase 
space of the discrete system. The function f which defines a discrete dynamical system 
is called a mapping, or a map, for short. The multipliers of the map f are the 
eigenvalues of the Jacobian of f. An area-preserving map is a map whose Jacobian 
has a unit determinant. 
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The set of iterates of f - x, /(x), f 2 (x), / 3 (x),..., f n (x) - as n becomes large is 
called the orbit of x relative to f; it captures the history of x as f is iterated. 

Two types of point have the simplest histories — fixed point, and periodic point. 
The point x is a fixed point of f if f(x) = x. A fixed point x is called stable, or 
elliptic, if all the multipliers of f at x lie on the unit circle; it is called unstable, or 
hyperbolic, otherwise. The point x is a periodic point of period n if f n (x ) = x. 
The least positive n for which f n (x) = x is called the period of x. The set of all 
iterates of a periodic point forms a periodic orbit. 


3.2 Analytical Techniques 

The first step in analyzing any dynamical system is to find its equilibria and their 
stability. For a second order discrete system, this amounts to solving a pair of non¬ 
linear algebraic equations. To take an example, the fixed points of quadratic map are 
the solutions to the following pair of algebraic equations: 

x = x cos a — (y — x 2 ) sin a 

y = x since + (y — a; 2 ) cos a: (2) 

There are two solutions: x = y = 0 and x = 2 tan(f), y = 2 tan 2 (f). 

Once the fixed points are found, the trace of the Jacobian of the map at each fixed 
point is computed to determine its linear stability. The stability condition is that the 
absolute value of the trace of the Jacobian is less than or equal to two [7]. For the 
quadratic map, the Jacobian is given by: 

( cos at 4- 2x sin a — sin a \ 
sin a — 2x cos a cos a J 

At (0,0), the trace of the Jacobian is 2 cos a. Since | 2 cos a | < 2, (0,0) is always 
stable. Similar stability analysis shows the other fixed point is always unstable. 


3.3 Qualitative Dynamics and their Geometry 

3.3.1 Types of Orbit 

In this section, I give a brief outline of the most important types of orbits in area¬ 
preserving maps. One important physical situation that gives rise to area-preserving 
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map is conservative systems of two degrees of freedom. Restricted to constant en¬ 
ergy surface, the 4-dimensional phase space becomes 3-dimensional. By the Liouville 
Theorem on integrable systems [1, Chapter 10], this 3-dimensional subspace is dif- 
feomorphic to a family of concentric 2-dimensional tori lying inside one another (see 
figure 4). A torus is called invariant if an orbit starting at a point on the surface of 
the torus stays on it forever. 



Figure 4: 3-dimensional phase space filled with 2-dimensional concentric tori. A 
cross-section E is drawn. 

Qualitatively, three special types of motion are important: (1) periodic motion 
lying on an invariant torus, (2) quasiperiodic motion that densely covers an in¬ 
variant torus, and (3) chaotic motion that wanders in a 3-dimensional volume of 
the phase space. 

We construct an area-preserving map by taking a 2-dimensional cross-section E of 
the family of tori. An orbit beginning in E returns to it after making a circuit around 
the torus. The successive intersections of the orbit with E defines a map from E onto 
itself. The map is area-preserving because the 3-dimensional flow is conservative [2, 
Appendix 31]. 

Now, how do the periodic and quasiperiodic orbits appear on E? When the motion 
is periodic, the orbit closes after a finite number of intersections with E, resulting in 
a finite number of points on E. If the motion is quasiperiodic, it will intersect E 
an infinite number of times, covering a closed curve densely. Finally, if the motion 
is chaotic, the orbit will intersect E erratically at many different places, eventually 
covering a finite area. 

Including the case of unbounded motion for the sake of completeness, we have the 
following correspondence between special types of orbit in the 3-dimensional phase 
space and their geometric manifestations on the corresponding 2-dimensional phase 
space of the map: 
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periodic motion 
quasiperiodic motion 
chaotic motion 
unbounded motion 


periodic points 
<=> invariant curve 
<=$■ chaotic region 
■$=> escape points 


3.3.2 Bifurcations: Qualitative changes in the phase portrait 

Two phase portraits are qualitatively equivalent if there exists a homeomorphism 
between them which preserves fixed points, periodic points, invariant curves, and 
their stability. Bifurcation is said to occur when the dynamical system goes through 
a qualitative change in its phase portrait as the control parameter is varied. I will 
focus on one important type of bifurcation: appearance and disappearance of periodic 
orbits. 

I begin with [12] which gives a complete classification of the generic bifurcations of 
periodic points for one-parameter area-preserving maps. Only bifurcations which are 
generic are classified because the pathological cases arise rarely, and are much harder 
to deal with. Let A and ^ be the multipliers of a periodic point of an area-preserving 
map. Meyer has shown that generic bifurcations occur when A = n th root of unity 
where n = 1, 2, 3, 4, and > 5. Let us consider each of these cases (see figure 5): 


Case 1 . A = 1 Extremal Bifurcation 


The (discrete) flow line forms a cusp at the bifurcation point and as the multiplier 
passes through +1, there is a sudden birth of a pair of elliptic and hyperbolic closed 
orbits. Thus, a pair of fixed points are created for the map. This is also known as 
creation , saddle-node bifurcation, or tangent bifurcation. 


Case 2. A = —1 Transitional Bifurcation 


As the multiplier passes through -1, the stable elliptic periodic point loses its 
stability, i.e., turns hyperbolic, and a stable cycle of twice the period appears. Other 
names for this phenomenon are: period doubling, flip bifurcation, subtle division. 


Case 3. A = e 3 Phantom 3-kiss 
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The region of stability of the elliptic fixed point shrinks to zero as the hyperbolic 
points of an unstable period-3 cycle “kiss” at the origin. After the “kiss”, the fixed 
point turns elliptic again, and a new unstable period-3 cycle is emitted. Note the 
change in orientation of the triangular region around the elliptic point. The phantom 
3-kiss is often preceded by extremal bifurcations in a region a bit further away from 
the original elliptic fixed point, resulting in the formation of a pair of elliptic and 
hyperbolic period-3 points. 


Case 4. A = e « Phantom 4-kiss 


When the multiplier A is a 4-th root of unity, it can undergo two types of bifurca¬ 
tion. The first case, phantom 4-kiss, is similar to that of case 3 except that we now 
have hyperbolic points of a period-4 orbit “kissing”. Tangent bifurcation also occurs 
at a distance from the center elliptic fixed point, resulting in a an unstable period-4 
cycle. The second case belongs to the fifth type to be described next. 


Case 5. A = e p where p > 5 Emission 


As the multiplier passes through a p-th root of unity where p > 5, there is an 
emission of elliptic and hyperbolic sub-harmonics of k of the orbital frequency of the 
original elliptic orbit (shown in the figure for the case p = 5). The elliptic points of 
the stable period-p cycle alternate with the hyperbolic points of the unstable one near 
the vertices of a regular p-gon with center at origin. Thus, emission is responsible for 
the creation of “island chains” in the phase portraits of Figure 1. Meyer called this 
type of bifurcation the generic p-bifurcation. 


4 APPROACH TO THE CONTROL PROBLEM 

4.1 How to start the numerical experiment? 

Elliptic fixed points are good places to start. We expect that the orbits near an 
elliptic fixed point, where the linear terms of the map dominate, will be mostly 
invariant curves. We then search radially outward until we encounter island chains, 
and eventually chaotic regions. 
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Figure 5: Five Generic types of Bifurcation Geometry. A periodic point bifurcates 
whenever its multipliers pass through an n-th root of unity. The bifurcation diagram 
shows the loci of fixed points as the parameter t is varied. 

4.2 How to decide what experiment to try next? 

Knowing the generic bifurcation patterns is valuable for controlling numerical exper¬ 
iments. To begin with, it is difficult to locate the value of the control parameter 
at which bifurcation occurs: it is of probability almost zero that a randomly chosen 
point in the (x, y , a) space will be the bifurcation point. But the pattern of flow near 
a periodic point just before and after the bifurcation occurs in a finite range of the 
control parameter; hence the pattern is easier to detect during the experiments. 

Once a given flow pattern is found to match some parts in our library of bifurcation 
geometries, it will give us strong evidence that the corresponding bifurcation exists, 
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and we should be able to locate the rest of the flow patterns as given by the generic 
bifurcation. The pre-stored knowledge about these bifurcations gives us the complete 
information about what geometric objects, and approximately where in the control 
parameter space to look for. 

To take an example, consider the phantom 3-kiss seen in figure Id. The local flow 
pattern around the fixed point matches that in figure 5. According to the bifurcation 
pattern, the regular region around the stable fixed point will shrink in size, becoming 
an unstable fixed point; eventually, a new stable fixed point is born. So, we should 
expect to see figure If at some a slightly greater than two. 


4.3 How to decide when to terminate the experiments? 

Besides imposing a strong constraint on what can be expected to happen in the phase 
portrait, the generic bifurcations also provide an answer to the problem of termina¬ 
tion: a simulation experiment is incomplete unless all the major qualitative features 
in the phase portrait can be explained by this finite list of local generic bifurcations. 
An example is the change of stability of a fixed point. Suppose we have numerically 
located the 3-island chain and the center elliptic point at some a = a 0 as in figure 
Id. We know that the family of phase portraits is yet incomplete because we expect 
a phantom 3-kiss bifurcation to occur. In particular, we need to try at least two more 
experiments to obtain two phase portraits: first, at a = otj when the triangular region 
is flipped, and second, at a 2 € (cco, c*i) when the region becomes vanishingly small, 
indicating instability of the fixed point. 


5 SOLVING THE INTERPRETATION PROBLEM 


The Interpretation Problem consists of the following sub-problems: 

1. Orbit Type. How can one recognize the orbit type - a O-dimensional finite point 
set whose elements are encountered repeatedly, a 1-dimensional smooth curve, 
or a 2-dimensional region - of a set of iterates? 

2. Clustering. How can one determine the number of islands in an island chain? 
This number gives the period of the enclosed periodic point. 
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3. Area and Centroid. How can one estimate the centroid and area enclosed by 
the curve? The centroid is a good approximation of the location of the enclosed 
periodic point. The area gives a measure of saliency of the island chain. 

4. Shape. How can one recognize the shape of a curve? For example, is it a 3-sided 
figure resembling a triangle? 


In the following, I will show how these four problems can be solved by applying 
techniques from computational geometry and computer vision. 


5.1 Euclidean Minimal Spanning Tree 

Given n points in the Euclidean plane, a Euclidean minimal spanning tree (EMST) is 
a tree that interconnects all the points with minimal total edge length. It turns out 
that EMST provides us with the fundamental data structure to solve the first three 
sub-problems. 

Figure 6 shows the steps of processing. First, the program computes a EMST from 
the input point set. The current implementation uses the Prim-Dijkstra algorithm 
which has a run time complexity of 0(n 2 ) where n is the number of points [3]. It takes 
about 100 seconds on the Symbolics to compute a EMST for a set of 512 points. A 
faster algorithm due to Shamos [13] which makes use of the Voronoi diagram is being 
implemented. The new algorithm can achieve the optimal lower bound of 0(n log n) 
time. Figure 7 shows the EMSTs constructed for three sample point sets. 

Second, the program detects clusters in the EMST by looking for edges in the tree 
that are significantly longer than nearby edges. Such edges are called inconsistent 
[17]. The following criterion suggested by Zahn is used to determine if an edge is 
inconsistent. An edge is inconsistent if (1) its length is more than two times the 
average of nearby edges, and (2) its length is more than two standard deviations 
larger than the average of the lengths of nearby edges (where the statistics are taken 
from the set of all nearby edges). A point P is nearby point Q if P is connected to 
Q in the EMST containing two or fewer edges. Inconsistent edges are then deleted. 
This breaks up the EMST into connected sub-components which are collected by a 
depth-first tree walk. Figure 8 shows the result of the clustering procedure on the 
input point set of figure 7b. The procedure finds three inconsistent edges. For the 
other two point sets, no inconsistent edge is detected. The output of the clustering 
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INPUT: A POINT SET 



Figure 6: Main processing steps for solving the interpretation problem 

procedure is a collection of connected sub-components each of which represents a 
sub-tree of the original EMST. 

Third, for each sub-tree of the EMST, the program examines the degree of each 
of its nodes, where the degree of a node is the number of nodes connected to it in 
the sub-tree. For a smooth curve, the EMST consists of two terminal nodes of degree 
one; the rest, degree two. For a point set that fills an area, its corresponding EMST 
consists of many nodes having degree three or higher. 

Finally, to compute the area and centroid of the region bounded by a curve, 
the program generates an ordered sequence of points from the EMST, and spline- 
interpolates the sequence to obtain a smooth curve. The smooth curve is encoded 
using chain coding [6]. Straightforward algorithms are then applied to compute the 
area and centroid. 
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Figure 8: Result of clustering. Three inconsistent edges are detected in the point set 
of figure 7(b). 

5.2 Scale Space Image 


The shape of a curve, as we have seen in section 4, is an important clue to determine 
what type of bifurcation is occurring. The number of sides, the locations and types 
of extremal curvature — these are two of the most important geometric properties to 
recognize. 

The recognition algorithm is based on the method of scale space image introduced 
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by Witkin [16]. First, a curve is parameterized by C(s) = (x(s),y(s)) where s is 
the arc length along the curve. The two functions x(s) and y(s) are periodic for a 
closed curve, and they can be computed from the chain code representation. Second, 
x(s) and y(s) are smoothed by the Gaussian and its first two derivatives at multiple 
spatial scales. The scales are separated by one octave: 4, 5, 7, 11, 15, and 22, which 
are results of successively multiplying 4 by y/2 and rounding. Third, the curvature 
function k(s) is computed by the formula: 

_ xy - yx 
(x 2 + 1 / 2 )! 

Finally, the zero-crossings of k(s), and the signs of k(s) are computed to determine 
the locations and type of the extrema. 



Figure 9: Type and locations of extrema of curvature. Spatial scales: 4, 5, 7, 11, 15, 
and 22. • = maximum curvature; x — minimum curvature. 

An example of the operation of the recognition algorithm is shown in Figure 9. 
The input is an orbit segment of 256 points which seems to fill a 3-sided curve. The 
positions of the extrema of curvature found at various smoothing scales are marked: 
• for maximum curvature; x for minimum curvature. 
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6 Summary 


In this paper, I have studied the task of qualitative analysis of nonlinear area¬ 
preserving map by numerical experiments. I have also described how to approach 
the two major problems in automating the experimenting process: (1) experiment 
control, and (2) result interpretation. The basic idea is that knowledge of qualitative 
dynamics and bifurcations provides a strong constraint on the type of behavior pos¬ 
sible. Finally, I have described a program which solves the interpretation problem by 
using techniques from computational geometry and computer vision. 
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